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1.  INTRODUCTION 


This  report  describes  the  work  performed  under  Contract  DAAH01-82-C-0519, 
"Tactical  Missile  Base  Flow,"  which  addresses  Task  II  of  Technical  Requirement 
0119,  of  October  21,  1980,  The  purpose  of  the  Task  II  effort  is  to  investi¬ 
gate  the  feasibility  of  applying  advanced  Navier-Stokes  numerical  simulation 
techniques  to  the  complex  flow  that  develops  near  the  base  of  a  jet-propelled 
tactical  missile  or  a  rocket-assisted  artillery  projectile.  This  complex 
flow  has  been  analyzed  in  the  past  with  the  help  of  ad-hoc  "component"  models, 
which  provide  relatively  low-cost,  approximate  solutions  for  engineering  pur¬ 
poses  [Anon.  (1969)].  The  rapid  gains  being  made  in  numerical  modeling  based 
on  the  Navier-Stokes  equations,  together  with  dropping  costs  of  computer 
utilization,  lead  naturally  to  attempting  the  application  of  these  techniques 
to  the  base  flow  problem.  Specifically,  successful  computations  of  base  flow 
with  simulated  jet  (modeled  by  a  solid  body  of  appropriate  shape)  made  at  the 
NASA  Ames  Research  Center  [Deiwert  (1981),  (1982)]  represent  a  leading  step  in 
this  direction,  and  open  the  way  to  numerical  simulation  of  the  fully  coupled 
near  wake  with  gas  jet. 

The  technical  challenge  is  posed  by  the  complexity  of  the  base  flow,  which 
involves  the  Interaction  between  an  underexpanded  propulsive  jet  and  the 
recirculating  flow  in  the  near  wake.  In  turn,  the  wake  flow  may  affect  the 
afterbody  flow,  causing  separation  of  the  boundary  layer  near  the  missile 
base,  depending  on  the  shape  of  the  body  and  the  base  pressure  that  develops 
as  a  result  of  the  jet-wake  interaction. x  For  almost  all  applications,  tur¬ 
bulent  flow  develops  in  the  shear  layers  that  separate  the  jet,  the  recir¬ 
culating  base  flow,  and  the  free  stream. 

The  approach  adopted  to  meet  this  challenge  consisted  of  restricting  the 
first  attempt  to  axisymmetric  flow  with  a  cold  propulsive  jet  of  the  same  gas 
as  the  free  stream  (both  gases  being  thermally  and  calorically  perfect) ,  and 
using  the  Lockheed  Viscous  Implicit  Solver  (LVIS),  a  Navier-Stokes  computer 
code  already  tested  in  a  number  of  other  applications  [see  Reklis,  et.  al. 
(1983)].  The  effort  concentrated  on  adding  to  LVIS  the  capability  to  treat 
the  base  region  (i.e.,  to  work  with  an  L-shaped  computational  domain),  and  to 
incorporate  two  different  turbulence  models,  of  the  so-called  "two-equation” 
type,  namely  the  k-epsilon  and  the  k-W  models  (see,  for  instance,  [Launder  and 
Spalding  (1972)]).  Validation  of  numerical  simulations  is  always  an  important 
aspect  of  their  application  to  practical  problems,  and  in  this  respect  the 
approach  adopted  was  to  establish  a  "sealed  envelope"  test,  according  to  which 
the  code  would  be  applied  in  a  true  predictive  fashion  to  wind-tunnel  tests 
performed  at  AEDC,  the  results  of  which  would  be  unavailable  to  the  contractor 
performing  the  numerical  simulation.  Thus,  the  LVIS  code  was  modified  as 
required,  and  applied  to  the  specific  AEDC  test  conditions,  using  both  the  k- 
epsilon  and  k-W  turbulence  models.  Comparison  of  these  results  with  the  AEDC 
experimental  data  will  be  performed  by  MICOM,  and  the  conclusions  from  this 
exercise  will  be  forthcoming. 

This  final  report  includes  the  equations  and  boundary  conditions  used  in 
the  computer  code  (Section  2),  the  numerical  solution  technique  (Section  3), 
and  a  description  of  the  results  of  the  numerical  simulation  computer  runs 
made  with  the  k-epsilon  and  k-W  turbulence  models  (Section  4).  The  detailed 
data  from  the  numerical  simulations  are  available  in  magnetic  tape.  The 
instructions  for  using  the  computer  code  are  given  in  a  Users  Manual  that  is 
contained  in  a  separate  volume. 
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2.  EQUATION  AND  BOUNDARY  CONDITIONS 

2.1  Governing  Equations 

2.1.1  Flow  Equations.  The  turbulent  flow  over  the  forebody  and  in  the 
base  region  is  assumed  to  be  governed  by  the  classical  Navier-Stokes  equations 
with  effective  values  of  viscosity  and  thermal  conductivity  that  are  the  sum 
of  laminar  and  turbulent  contributions 


(2.1a) 


(2.1b) 


where  p  denotes  the  viscosity,  K  the  thermal  conductivity,  Cp  the  heat  capac¬ 
ity,  Pr  the  Prandtl  number,  and  the  subscripts  1  and  t  refer  to  molecular 
(laminar)  and  turbulent  quantities,  respectively.  In  the  present  work,  the 
eddy  viscosity  Is  obtained  from  either  one  of  two  multi-equation  turbulence 
models,  and  the  turbulent  Prandtl  number  is  assumed  to  be  a  constant. 

The  numerical  solution  to  the  described  Navier-Stokes  equations  is  carried 
out  in  a  general  curvilinear  coordinate  system.  The  curvilinear  coordinate 
transformations  and  equations  can  be  found  in  [Pulliam  and  Steger  (1980)  or  in 
Thomas  and  Lombard  (1970)].  The  simplification  of  the  equations  for  axlsym- 
metric  flow  is  described  by  [Nietubicz,  et  al.,  1979].  In  the  present  work, 
the  full  set  of  viscous  terms  are  retained  (see,  for  example,  [Thomas  (July, 
1979)],  including  the  cross-derivative  terms;  whereas,  these  terms  were 
omitted  in  [Pulliam  and  Steger  (1980)]  and  in  [Nietubicz,  et  al.  (1979)]. 

2.1.2  Turbulence  Model  Equations.  In  high  Reynolds  number  flow  regions 
away  from  walls,  we  employ  a  generalized  two-equation  turbulence  model  that 
embodies  both  the  uncorrected  k-epsilon  model  [Jones  and  Launder  (1972), 
Launder  and  Spalding  (1972),  and  Mace,  et  al.,  1981]  and  the  k-W  model 
[Launder  and  Spalding  (1972)  and  Spalding  (1972)]. 

The  generalized  model  equations  in  dimensionless  form  are 


(pk)t  +  V* PVk  -  Re  V*(Dk7k)  - 

sk  -  P[^«  -  cDk  r®"t] 


(pr)t  +  7*pvr  -  Re  7*(Dp7 T)  -  Sr 


(2.2a) 


(2.2b) 


3r  -  p[cprkapr  Ar  $  +  c’rkapr  r0pr(Vco)2  -  cprkV  Ar] 


& 


where  k  denotes  the  turbulent  kinetic  energy  (T.K.E.),  T  represents  the  second 
turbulence  model  parameter  (e'psilon  or  W) ,  Dfc  and  Dp  are  diffusion  coeffi¬ 
cients  defined  in  terms  of  the  molecular  and  eddy  viscosities  end  "exchange 
coefficients"  a  that  are  empirical  constants 

Vt 

Dk  -  a •  +  —  (2.3a) 


(2.3b) 


and  where  p  is  the  density,  V  is  the  flow  velocity,  9  is  the  classical  dissi¬ 
pation  function,  which  can  be  represented  in  Cartesian  tensor  notation  as 


a)  =■  J \kV  |  is  the  magnitude  of  the  vorticity  vector,  vt  is  a  Reynolds  number- 
scaled  turbulent  kinematic  viscosity 

^  -  Cv  k^  (2.4) 

and  the  eddy  viscosity  is  given  by 


Mt  -  Re  p  vt  (2.5) 

All  quantities  in  the  above  equations  are  dimensionless,  with  velocity, 
distance,  density,  and  viscosity  normalized  by  some  reference  value  for  each, 
time  non-dimensionalized  by  the  ratio  of  reference  length  to  reference  veloc¬ 
ity,  turbulent  kinetic  energy  normalized  by  the  square  of  the  reference  veloc¬ 
ity,  T  non-dimensionalized  by  the  appropriate  combination  of  reference  length 
and  velocity,  and  the  Reynolds  number  Re  is  formed  from  the  reference  values 
of  density,  velocity,  length,  and  viscosity.  The  reference  conditions  used  in 
the  present  baseflow  analysis  are  the  freestream  speed  of  sound,  density, 
viscosity,  temperature,  and  pressure;  and  the  exit  radius  of  the  nozzle  (see 
also  Section  4). 

The  remaining  coefficients  and  exponents  in  the  foregoing  turbulence  model 
equations  are  constants.  Dimensional  analysis  determines  uniquely  the  values 
of  the  exponents  in  terms  of  the  dimensions  of  the  second  turbulence  parameter, 
T.  If  we  represent  the  latter’s  dimensions  in  units  of  length  L  and  the  time 
T  follows 


V*wv«r vmnrywvr  yT5r*3r*jr-j 


i*v  .4 


r  LaT.  U'*J 


[r]~L"T" 


then  the  appropriate  values  for  the  two  turbulence  models  are 
k-epsilon  Model:  n  ■  2,  m  ■  -3 


(2.6) 


k-W  Model: 


n  *  0,  m  ”  -2 


and  the  various  exponents  a,  p  in  Eq*s  (2.2)-(2.5)  are  given  by 

_  _  n+2m  _  1 

v  2(n-hn)  v  n-hn 


^  3n+2m 
2  (n+m) 


n+m 


n 

2 (n+m) 


n+2m 
2  (n+m) 


2 (n+m) 


m  I  —  _ 

•  T  n+m 


(2.7a) 


(2.7b) 


(2.7c) 


(2.7d) 


(2.7e) 


Table  2.1  gives  the  values  of  the  remaining  constants  in  Eq's  (2.2 )—( 2 . 5 ) 
for  which  those  equations  reproduce  either  the  k-epsilon  or  the  k-W  turbulence 
model  equations.  For  reference,  the  table  also  Indicates  the  equivalent  nota¬ 
tion  that  normally  has  been  used  for  the  coefficients  in  previous  work  dealing 
with  either  turbulence  model.  For  example,  the  symbol  C\  previously  has  been 
used  in  the  k-W  model  literature  to  denote  C'pp,  the  empirical  coefficient 
of  (Va))2  in  the  second  production  term  on  the  R.H.S.  of  Eq.  (2.2b). 

The  turbulence  model  equations  (2.2)-(2.5)  apply  only  in  high  Reynolds 
number  flow  regions  far  from  walls.  In  low-speed  flow  regions,  such  as  within 
the  boundary  layers  on  the  missile  side  wall  and  on  the  solid  part  of  the  base, 
one  either  must  add  further  "wall  correction"  terms  to  the  model  equations,  or 
must  employ  some  other  near-wall  turbulence  model.  Several  investigators  have 
introduced  "wall  correction”  terms  into  the  k-epsilon  model  equations  and 
applied  the  expanded  models  successfully  to  simple  boundary  layer  flows  [Jones 
and  Launder  (1972),  Chien  (1982)].  However,  the  behavior  of  such  correction 
terms  has  yet  to  be  investigated  for  complex  flows  with  vortices,  embedded 
shock  waves,  etc.  Therefore,  we  have  elected  to  follow  the  approach  of 
[Arora,  et  al.,  (1982)],  who  used  an  algebraic  mixing  length  model  for  the 
eddy  viscosity  1  the  near-wall  region  of  the  boundary  layers,  and  who  applied 
the  turbulence  model  equations  (2.2)-(2.5)  only  outside  this  region. 


Table  2.1 


Constants  In  Generalized  Turbulence  Model  Equations  (2.1) 


Constant 


k-epsllon  Model 
[Mace,  et  al  (1981)] 


k-W  Model 
[Spalding  (1972)] 


Value  of 
Constant 


Corresponding 

Symbol 


Value  of 
Constant 


Corresponding 

Symbol 


For  smooth  walls  [Arora,  et  al.  (1982)]  use  the  van  Driest  eddy  viscosity 
formulation 


vt  *  (0.4  Dy)2  $ 
D  ■  1  -  exp(-E) 


(2.8) 


where  y  is  the  distance  normal  to  the  wall.  For  generality,  we  have  employed 

1/2 

the  square  root  of  the  dissipation  function  <*  instead  of  the  wall-normal 
gradient  of  tangential  velocity  as  used  in  [Arora,  et  al.  (1982)];  the  two 
become  equal  in  the  neighborhood  of  a  wall. 

In  [Arora,  et  al.  (1982)],  the  exponent  E  of  the  damping  factor  D  is  given 


A+  %/  A+ 


t  Pw  Hw^jr 


my/2 


(2.9) 


where  y+  is  the  distance  from  the  wall,  normalized  by  the  viscous  length 
scale,  which  is  defined  as  the  ratio  of  the  kinematic  viscosity  at  the  wall  to 
the  friction  velocity;  x  is  the  local  shear  stress,  and  A+  is  the  damping 
constant. 

An  alternative  expression  is  the  formula  used  for  the  inner  near-wall 
region  in  the  Baldwin-Lomax  algebraic  turbulence  model  [Baldwin  and  Lomax 
(1978)] 


r  *1/2i 

p  <& 

w  w 


(2.10) 


In  the  present  work,  as  in  [Arora,  et  al.,(1982)],  the  edge  of  the  inner 
region  is  to  be  positioned  just  outside  the  laminar  sublayer,  y+  15. 

Boundary  conditions  for  the  two  turbulence  model  parameters  at  this  edge  are 
obtained  from  the  dual  condition  that  the  eddy  viscosity  must  be  continuous 
across  the  edge,  and  that  the  turbulent  kinetic  energy  k  is  assumed  to  be  in 
local  equilibrium.  The  continuity  condition  gives  one  equation  between  k  and 
Tat  the  inner  region  edge  by  equating  the  two  expressions  (2.4)  and  (2.8)  for 
vx  and  the  local  equilibrium  condition  yields  a  second  equation  by  setting 
equal  to  zero  the  source  term  ^k  on  the  R.H.S.  of  the  T.K.E.  Equation  (2.2a). 

A  general  three-dimensional  coordinate  transformation  is  applied  to  the 
unsteady  Navier-Stokes  equations  [Thomas  and  Lombard  (1979)]  and  the  trans¬ 
formed  equations  are  simplified  for  axisysmmetric  flow  [Nietubicz,  et  al., 
(1979)].  The  same  transformations  and  simplifications  are  applied  to  the  tur¬ 
bulence  model  Equations  (2.2)  -  (2.10).  If  we  represent  the  transformation  in 
the  form 


i 

I 


x,  y,  z,  t- 


^  £»  h*  C»  t 


with  t  ■  t,  and  take  t|  as  the  meridional  angle  in  axisymmetric  geometry,  then 
the  transformed  turbulence  model  equations  can  be  written  in  vector  form  as 


(jpQ)T+  [(Jpv-  v^)q]  j.  +  [(jpv.  7c)q]c  "  r®1[di^f,  +  °3 Pr] 
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and  where  J  -  3(x,y ,z) /  d(  n,  O  is  the  Jacobian  of  the  inverse  transformation, 
and  D  is  a  diagonal  diffusion  coefficient  matrix 
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3.  NUMERICAL  SOLUTION  TECHNIQUE 


The  numerical  solution  to  the  described  equations  is  carried  out  using  the 
Beam-Warming  time-linearized  implicit  approximate-factorization  (ADI)  scheme 
[Beam  and  Warming,  (1978)].  When  this  scheme  is  used  with  algebraic  turbu¬ 
lence  models,  the  usual  practice  is  to  treat  both  laminar  and  turbulent  trans¬ 
port  properties  as  locally  constant  over  a  time  step  At,  and  then  to  update 
these  quantities  at  the  end  of  the  step.  Details  of  the  implicit  scheme 
applied  in  this  fashion  in  curvilinear  coordinates  are  given  in  [Pulliam  and 
Steger,  (1980)].  For  multi-equation  turbulence  models  such  as  employed  here, 
one  could  solve  the  entire  system  of  flow  and  turbulence  model  equations 
directly  by  the  implicit  numerical  scheme,  but  this  would  be  both  cumbersome 
and  computationally  inefficient.  Instead,  we  invoke  a  few  simple  approxima¬ 
tions  that  do  not  affect  the  accuracy  of  the  steady  state  solutions  of  inter¬ 
est  here,  and  that  both  simplify  the  scheme  and  improve  the  computational 
efficiency.  These  approximations  and  the  resulting  computational  procedure 
are  described  below. 


3.1  Flowfield  Difference  Equations.  The  flowfield  equations  depend  on  the 
turbulence  model  parameters  k,  I*  only  through  the  eddy  viscosity  p^..  As  with 
algebraic  turbulence  models,  we  neglect  this  dependence  to  first  order  by 
treating  ^  as  locally  constant  over  a  step  At.  This  allows  the  flowfield  and 
turbulence  model  equations  to  be  solved  in  tandem,  rather  than  simultaneously. 
The  flow  equations  are  solved  first  in  the  usual  way  [Thomas  and  Lombardy 
(1979),  Pulliam  and  Steger  (1980)]  to  obtain  the  time  differences  AP,  AV  and 
the  updated  flow  variables  Pn+*,  \fc+l,  etc.  These  values  then  are  used  in  the 
subsequent  solution  of  the  turbulence  model  equations  by  essentially  the  same 
implicit  scheme. 


3.2  Turbulence  Model  Difference  Equations.  We  employ  the  first-order 
Euler-implicit  time  differencing  scheme.  The  time- linearized  difference 
equations  for  Eq.  (2.11)  then  can  be  written  In  the  form 
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in  which  I  is  the  identity  operator  (matrix)  and  £Q  is  the  only  unknown. 

Here,  we  have  employed  classical  difference  operator  notation  A,  6,  p6  with  a 
subscript  on  each  spatial  operator  to  distinguish  the  direction  in  which  that 
operator  acts.  For  example,  for  a  mesh  function  f”  ,  the  spatial  operators 
are  defined  as  follows  J 


Lf  -  f 


1+1, i  "  fi,j 


»*Af  “  (fi+l,j  "  fi"l,j)/2 

n  n 

*  fi+l/2  “  fi-l/2 


(3.2) 


and  A  without  a  subscript  represents  the  forward  time  difference  operator 


(3.3) 


Following  [Beam  and  Warming  (1978)],  the  diffusive  cross-derivative  terms 
have  been  treated  in  explicit  fashion  in  Equation.  (3.1).  For  later  reference, 
we  note  that  difference  equation  (3.1a)  does  not  contain  artificial  smoothing 
terms. 

To  simplify  the  solution  of  the  difference  equations  (3.1),  we  perform  an 
approximate  factorization  of  the  implicit  L.H.S.  operator  into  three  parts. 

Two  of  these  are  the  usual  one-dimensional  factors  for  each  coordinate 
direction,  and  the  third  involves  only  the  source-term  Jacobian  matrix  dS/dQ. 
The  factored-operator  form  of  Equation  (3.1)  is 
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The  implicit  operator  inversion  sequence  thus  Involves  two  successive 
block-tridiagonal  matrix  inversions,  one  for  each  spatial  direction,  followed 
by  a  block-explicit  matrix  inversion  locally  at  each  grid  point,  where  each 
block  has  only  a  dimension  of  2.  In  each  of  the  spatial  block-tridiagonals, 
the  two  equations  (for  the  two  components  of  Q)  are  coupled  only  through  the 
turbulent  viscosity  Pj..  The  block-tridiagonal  then  is  reduced  to  a  set  of  two 
easily-inverted  scaler  tridiagonal  equations  upon  approximating  as  locally 
constant  over  a  time  step,  which  is  the  same  approximation  employed  in  treating 
the  flow  equations  themselves.  Note^that  this  approximation  also  renders  the 
diffusion  terms  T^  and  T2  linear  in  Q,  and  their  Jacobian  matrices  take  the 
simple  form 
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Numerical  experiments  have  revealed  that  the  described  algorithm  for  the 
turbulence  model  equations  is  both  stable  and  accurate  in  regions  where  the 
turbulence  level  is  significant,  such  as  in  a  boundary  layer,  shear  layer,  or 
jet  mixing  layer.  However,  numerical  difficulties,  usually  in  the  form  of 
slightly  negative  values  of  the  turbulence  energy  k,  often  arise  In  low- 
turbulence  regions  where  pt  in  Equation  (2.1)  is  insignificant.  An  example  is 
the  approximately  uniform  near-freestream  flow  region  radially  far  from  the 
missile  afterbody  wall  in  Figure  3.1,  where  k  and  p^  are  very  small.  These 
difficulties  have  been  overcome  by  employing  upwind^rather  than  central  dif¬ 
ference  operators  In  the  advectlve  terms  involving  V  on  both  right-hand  and 
left-hand  sides  of  the  difference  equation  (3.4).  The  Implementation  for  the 
turbulence  model  equations  follows  that  described  in  [Reklis  and  Thomas 
(1981 )J  for  the  Euler  equations. 

3.3  Boundary  Conditions.  Boundary  conditions  for  the  flow  variables  and 
turbulence  model  parameters  at  the  various  boundaries  depicted  in  Figure  3.1 
were  determined  and  applied  as  described  below. 

Implicit  boundary  conditions  at  the  outflow  plane  and  at  the  symmetry  axis 
were  computed  as  described  In  [Thomas  and  Lombard  (1979)].  The  remaining 
boundary  conditions  were  applied  in  explicit  fashion;  that  is,  they  were  main¬ 
tained  as  locally  constant  over  a  time  step  and  then  updated  at  the  end  of  the 
step.  No-slip  adiabatic  wall  boundary  conditions  were  used  at  the  missile 
sidewall  and  along  the  solid  part  of  the  base.  The  pressure  and  temperature 
at  each  grid  point  along  the  wall  were  updated  by  setting  them  equal  to  the 
values  at  the  adjacent  interior  point.  Boundary  conditions  for  the  two  tur¬ 
bulence  model  parameters  are  not  applied  at  the  wall,  but  at  the  edge  of  the 
near-wall  region  within  which  the  algebraic  eddy  viscosity  model  (2.8)  is 
used.  The  edge  of  the  near-wall  (y+«15).  Boundary  values  of  k  and  r  at  the 
edge  were  obtained  from  the  dual  requirement  that  p^  be  continuous  and  that  k 
be  in  local  equilibrium,  as  explained  earlier. 

Boundary  conditions  at  the  inflow  boundary  and  shock  were  derived  from 
independent  numerical  solutions  for  the  external  flow  over  the  missile 


forebody  and  for  the  internal  nozzle  flow,  and  were  held  fixed  during  the 
baseflow  computation.  The  forebody  flow  solution  was  computed  with  a  2-D  axi- 
symmetric  flow  version  of  the  LVIS  (Lockheed  Viscous  Implicit  Solver)  program 
described  in  [Reklis,  Conti,  and  Thomas  (1983)].  The  shape  of  the  outer  bow 
shock  wave  and  the  flow  conditions  at  the  shock  were  obtained  from  the  same 
solution  for  a  fictitiously  lengthened  body  extending  to  the  outflow  plane  of 
Figure  3.1. 

Inflow  boundary  conditions  at  the  nozzle  exit  were  obtained  from  a 
Navier-Stokes  internal  flow  solution  using  the  N0ZL3D  code  [Thomas  (Sept., 

1979)  and  (1980)  and  Thomas  and  Neier  (1980)]  in  a  recent  version  that  includes 
an  axisymmetric  flow  option.  The  geometry  of  the  nozzle  is  shown  in  Figure 
3.2,  along  with  the  curvilinear  grid  used  from  the  N0ZL3D  code  solution. 

Plots  of  selected  exit-plane  flow  profiles  from  the  latter  solution  are 
displayed  in  Figures  3.3  to  3.5  to  demonstrate  the  significant  radial  nonuni¬ 
formity  of  the  jet  produced  by  this  particular  nozzle.  For  example,  the  exit- 
plane  pressure  is  about  20  percent  higher  at  the  axis  of  symmetry  than  at  the 
nozzle  lip.  In  these  figures,  the  radius  z  is  normalized  by  the  nozzle  exit 
radius;  pressure  and  temperature  are  normalized  by  the  corresponding  stagna¬ 
tion  chamber  conditions,  and  velocities  by  the  stagnation  sound  speed. 

The  described  numerical  solutions  for  the  external  forebody  flow  and  the 
Internal  nozzle  flow  employed  the  Baldwin-Lomax  turbulence  model  [Baldwin  and 
Lomax,  (1978)].  Because  this  model  is  algebraic,  it  cannot  directly  yield 
Inflow  boundary  conditions  on  the  turbulence  model  parameters  k  and  T.  It  was 
necessary  to  derive  the  latter  boundary  conditions  with  the  same  technique 
described  earlier  for  determining  the  values  of  k  and  T  at  the  edge  of  the 
near-wall  region  of  the  boundary  layer.  In  the  case  of  the  inflow  boundary, 
the  known  value  of  from  the  Baldwin-Lomax  model  and  the  assumption  that  k 
be  in  local  equilibrium  at  each  point  allow  deduction  of  local  values  for  k 
and  T.  This  procedure  made  it  necessary  to  employ  Equation  (2.10)  for  the 
damping  factor  exponent,  which  is  consistent  with  the  Baldwin-Lomax  model,  in 
order  to  ensure  continuity  of  the  k  and  T  distributions  in  the  neighborhood  of 
the  inflow  boundaries. 

The  same  approach  was  used  to  derive  boundary  conditions  for  k  and  r  at 
grid  points  just  inside  the  bow  shock  boundary. 
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4.  RESULTS  OF  TEST  CASES 


4.1  General.  The  test  cases  computed  consisted  of  verification  of  the 
turbulence  models  in  a  boundary  layer,  and  simulation  of  the  flowfield 
corresponding  to  the  AEDC  wind-tunnel  experiments.  In  order  to  make  this  a 
true  prediction,  the  results  of  the  AEDC  measurements  have  not  been  released 
as  of  this  writing;  comparison  of  the  present  results  with  that  data  will  be 
made  by  MICOM  in  the  future. 

Computations  were  made  on  a  VAX  11/780  minicomputer  and  CRAY  Is  mainframe. 
The  same  code  (with  minor  modifications)  was  used  in  both  computers,  and  on 
occasion  the  computation  was  moved  from  one  computer  to  the  other  for  opera¬ 
tional  convenience.  Many  of  the  runs  were  of  experimental  nature,  and  the 
cumulative  running  time  was  not  computed.  The  running  time  of  a  "production" 
run  can  be  estimated  by  noting  that  the  code  runs  approximately  17  millisec¬ 
onds  per  grid  point  per  time  iteration  on  the  VAX  computer. 

4.2  Check  of  Turbulence  Models.  The  correct  operation  of  the  turbulence 
models  was  checked  by  computing  the  flow  in  the  boundary  layer  on  a  cylinder 
aligned  with  the  free  stream.  This  case  was  selected  for  convenience,  since 
the  base  flow  code  is  set  up  to  compute  such  a  flow  on  the  aft  section  of  a 
missile.  Three  runs  were  made,  using  the  Baldwin-Lomax  algebraic  model,  the 
k-epsilon  model,  and  the  k-W  model,  respectively.  The  parameters  used  in  the 
two-equation  models  (see  Table  2.1)  were  those  corresponding  to  axisymmetric 
flow,  even  though  a  more  faithful  modeling  of  this  particular  flow  should  use 
the  parameters  recommended  for  two-dimensional  flows  with  planar  symmetry. 

The  axisymmetric  flow  parameters  were  used  because  in  the  computation  of  base 
flows  these  were  the  parameters  to  be  used  thorughout  the  flowfield,  including 
the  boundary  layer  on  the  side  walls  of  the  missile. 

All  computations  yielded  essentially  the  same  pressure  profile.  The 
longitudinal  velocity  profiles  are  shown  on  Figure  4.1.  The  k-W  and  algebraic 
models  produced  essentially  the  same  velocity  profile,  but  the  k-epsilon  model 
gave  a  fuller  velocity  profile  in  the  region  extending  from  about  2  percent  to 
15  percent  of  the  boundary-layer  thickness.  This  behavior  was  also  observed 
in  the  density  profiles,  shown  on  Figure  4.2.  The  turbulent  viscosity 
(defined  in  equation  2.1a)  is  shown  on  Figure  4.3.  Here  again  the  algebraic 
and  k-W  models  agree  well,  but  the  k-epsilon  model  shows  significantly  larger 
viscosity  in  the  region  extending  from  about  3  percent  to  42  percent  of  the 
bound ary- layer  thickness. 

The  reason  for  the  discrepancies  shown  by  the  k-epsilon  model  remains 
obscure,  but  the  general  results  of  these  tests  indicate  that  the  two-equation 
models  are  operational  and  yield  reasonable  solutions.  Further  judgment  on 
the  performance  of  these  models  is  beyond  the  scope  of  this  study. 
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4.3  AEOC  Test  Case 


4.3.1  Discretization  Grid.  The  grid  used  in  the  computations  is  shown  on 
Figure  4.4.  This  grid  was  generated  using  the  technique  described  in  [Thomas 
(Sept.  1979)],  which  is  based  on  the  solution  of  a  special  set  of  elliptic 
differential  equations.  A  computer  code  (RGRID)  that  implements  this  grid 
generation  technique  is  provided  separately  [Thomas  and  Neier  (1980)]. 

The  grid  used  for  the  test  cases  has  active  3216  nodes.  High  resolution 
is  provided  in  the  boundary  layer  on  the  side  and  base  of  the  missile,  and  in 
the  shear  layers  that  separate  the  jet  from  the  recirculating  flow  and  the 
free  stream.  In  order  to  economize  grid  points,  the  resolution  was  decreased 
in  the  upper  region  of  the  flowfield,  especially  above  z-8,  which  is  the  upper 
boundary  of  the  region  surveyed  in  the  AEDC  experiments.  Therefore,  care  must 
be  used  in  interpreting  the  computed  results  in  this  upper  region. 

4.3.2  Boundary  and  Initial  Conditions.  The  boundary  conditions  are 
described  in  detail  in  Paragraph  3.3.  The  inflow  boundary  conditions  on  the 
afterbody,  at  x— 16,  were  obtained  from  a  separate  forebody  solution,  and  they 
match  closely  the  boundary  data  from  the  AEDC  tests.  The  inflow  boundary  con¬ 
dition  at  the  nozzle  exit  was  obtained  from  a  nozzle  computation,  as  discussed 
in  Paragraph  3.3  and  displayed  in  Figures  3. 2-3. 5. 

The  initial  flowfield  to  start  the  computation  was  obtained  from  a  laminar 
flow  solution  computed  on  a  coarse,  rectangular  grid.  The  turbulence  parame¬ 
ters  were  estimated  roughly  by  interpolation  of  their  boundary  values. 

The  free-stream  parameters  used  in  both  the  k-epsilon  and  k-W  solutions 
were  as  follows: 

Mach  number-1.343 

Reynolds  number-92320 

Prandtl  number-0.7 

Ratio  of  specific  heats-1.4 

Static  temperature-265.3  K. 

The  values  of  the  constants  used  in  the  turbulence  models  are  given  in  Table 

2.1. 

4.3.3  Time-Asymptotic  Convergence.  The  code  was  run  with  the  k-W  model 
from  the  initial  conditions  to  a  dimensionless  time  (see  definition  in 
Paragraph  2.1.2)  of  approximately  28,  until  no  appreciable  changes  occurred  in 
the  flow  variables  throughout  the  flowfield.  Then  the  k-epsilon  model  was 
run,  using  as  initial  conditions  the  k-W  solution  as  a  dimensionless  time  of 
approximately  8.  The  k-epsilon  solution  was  run  to  a  dimensionless  time  of 
approximately  23,  which  is  roughly  equivalent  to  the  age  of  the  k-W  solution, 
because  of  its  later  starting  point.  This  solution  was  also  converged,  in  the 
sense  that  no  further  changes  were  observed  in  the  flow  variables. 
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4.3.4  Computed  Flowfield.  The  k-W  solution  will  be  discussed  first.  The 
flow  variables  are  nondimensionalized  as  described  in  Paragraph  -  2.1.2;  in 
particular,  the  coordinates  are  referenced  to  the  nozzle  exit  radius.  Figure 
4.5  shows  the  general  structure  of  the  flowfield,  as  indicated  by  velocity 
vector  plots.  The  nozzle  exit  extends  from  the  centerline  (z»0)  to  z*l;  the 
base  of  the  body  extends  to  z«5.  The  near  wake  exhibits  a  recirculation 
region,  which  extends  longitudinally  to  approximately  x*6.  Within  it  there  is 
a  dividing  streamline  that  separates  the  upper  and  lower  recirculating  flows, 
with  a  stagnation  point  on  the  base  of  the  missile  at  approximately  z*2.9. 

The  shear  layers  that  separate  the  slow  wake  flow  from  the  jet  and  the  exter¬ 
nal  stream  are  clearly  visible  in  Figures  4.5(a),  (b),  and  (c),  which 
Illustrate  their  evolution  and  merging  along  the  wake. 

Two  minor  anomalies  remain  unexplained:  the  notch  in  the  velocity  profile 
on  the  jet  side  of  the  shear  layer,  particularly  visible  in  Figure  4.5(a);  and 
the  velocity  defect  on  the  axis  of  symmetry,  which  may  be  related  to  the 
axial-symmetry  boundary  condition. 

Dimensionless  pressure  and  density  contours  are  shown  in  Figures  4.6  and 
4.7,  respectively.  The  density  contours  demarcate  clearly  the  “edge"  of  the 
jet,  as  Figure  4.7  shows. 

The  k-epsilon  solution  is  very  similar  to  the  k-W  one.  Figures  4.8(a), 
(b),  and  (c)  show  the  velocity  vector  plots,  and  Figures  4.9  and  4.10  show  the 
pressure  and  density  contours,  respectively.  One  difference  with  the  k-W 
solution  is  that  the  k-epsilon  recirculating  region  is  shorter  by  about  one 
nozzle  exit  radius.  Other  differences  are  the  greater  width  of  the  k-epsilon 
shear  layers,  and  the  more  nearly  linear  velocity  profile  of  its  jet,  particu¬ 
larly  downstream.  These  differences  are  consistent  with  the  larger  eddy  vis¬ 
cosity  predicted  by  the  k-epsilon  model,  as  illustrated  by  Figure  4.11,  which 
shows  profiles  of  the  turbulent  viscosity  (mu  sub  tee)  at  two  axial  stations. 

The  computed  results  discussed  in  this  section  are  available  on  magnetic 
tape,  both  on  the  original  grid  described  in  Paragraph  4.3.1  and  on  the  grid 
surveyed  in  the  AEDC  wind-tunnel  tests. 
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Figure  4 . I 
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Comparison  of  computed  boundary  layer  velocity  profiles  for 
three  turbulence  models:  Q)  Baldwin-Lomax  algebraic  model 
(2)  k-epsilon  model,  (3)  k-W  model. 
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